Exploration and Exploitation in Parkinson’s Disease: Computational Analyses

Authors

Affiliations

Björn Meder

Health and Medical University, Potsdam, Germany

Martha Sterf

Medical School Berlin, Berlin, Germany

Charley M. Wu

University of Tübingen, Tübingen, Germany

Matthias Guggenmos

Health and Medical University, Potsdam, Germany

Table of contents

• 1 Preamble
• 2 Load data
• 3 Computational Analyses
  • 3.1 Gaussian Process Upper Confidence Bound (GP-UCB) Model
    • 3.1.1 Learning Model
    • 3.1.2 Sampling Strategy
    • 3.1.3 Choice rule
    • 3.1.4 Model parameters
• 4 Model comparison
  • 4.1 Model comparison: Control
  • 4.2 Model comparison: PD+
  • 4.3 Model comparison: PD-
  • 4.4 Model-based classification of participants
• 5 Analysis of parameter estimates
  • 5.0.1 Generalization \(\lambda\)
  • 5.0.2 Exploration bonus \(\beta\)
  • 5.0.3 Random exploration \(\tau\)
• 6 Relations of model parameters to performance
  • 6.1 Generalization \(\lambda\)
  • 6.2 Exploration bonus \(\beta\)
  • 6.3 Random exploration \(\tau\)
  • 6.4 Regression: Perfomance and model parameters
  • 6.5 Model simulations
• 7 Supplementary Information
  • 7.1 Statistical analyses
  • 7.2 Supplementary computational results
• 8 Article figure
  • 8.1 Figure 3. Computational results
• 9 Relations of model parameters and clinical indicators
  • 9.1 Bivariate correlations
  • 9.2 Regression analyses
• 10 Model params of participants best explained by GP-UCB model
  • 10.0.1 Generalization \(\lambda\)
  • 10.0.2 Exploration bonus \(\beta\)
  • 10.0.3 Random exploration \(\tau\)
• 11 Session Information

# Housekeeping: Load packages and helper functions
# Housekeeping
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(message = FALSE)
knitr::opts_chunk$set(warning = FALSE)
knitr::opts_chunk$set(fig.align='center')

options(knitr.kable.NA = '')

packages <- c('gridExtra', 'BayesFactor', 'tidyverse', "RColorBrewer", "lme4", "sjPlot", "lsr", "brms", "kableExtra", "afex", "emmeans", "viridis", "ggpubr", "hms", "scales", "cowplot", "waffle", "ggthemes", "parameters", "rstatix", "magick", "grid", "cetcolor", "ggcorrplot")

installed <- packages %in% rownames(installed.packages())
if (any(!installed)) {
  install.packages(packages[!installed])
}

# Load all packages
lapply(packages, require, character.only = TRUE)

set.seed(0815)

# file with various statistical functions, among other things it provides tests for Bayes Factors (BFs)
source('statisticalTests.R')

# Wrapper for brm models such that it saves the full model the first time it is run, otherwise it loads it from disk
run_model <- function(expr, modelName, path='brm', reuse = TRUE) {
  path <- paste0(path,'/', modelName, ".brm")
  if (reuse) {
    fit <- suppressWarnings(try(readRDS(path), silent = TRUE))
  }
  if (is(fit, "try-error")) {
    fit <- eval(expr)
    saveRDS(fit, file = path)
  }
  fit
}


# Setting some plotting params
w_box          <- 0.2      # width of boxplot, also used for jittering points and lines    
line_jitter    <- w_box / 2
xAnnotate      <- -0.3

# jitter params
jit_height  <- 0.01
jit_width   <- 0.05
jit_alpha   <- 0.6

# colors for age groups
groupcolors    <- c("#d95f02", "#1b9e77", "#7570b3")
choice3_colors <- c("#e7298a", "#66a61e", "#e6ab02")

1 Preamble

This document provides R code for the statistical analyses and plots of the behavioral data reported in the article

Meder, B., Sterf, M. Wu, C.M, & Guggenmos, M. (2025). Uncertainty-directed and random exploration in Parkinson’s disease. PsyArXiv

All analyses are fully reproducible, with the R code shown alongside the results, and random seeds set to ensure identical outputs across runs. Full session info is provided at the end of the document. All materials, including this document and all data, are available at:

https://github.com/charleywu/gridsearch_parkinsons

2 Load data

There are two files with behavioral data: data_gridsearch_Parkinson.csv

The behavioral data are stored in

• data_gridsearch_parkinson.csv, which contains the behavioral data from rounds 1-9 from the task
• data_gridsearch_subjects.csv, which contains participant information.

These files are combined to data frame dat, which includes the following variables:

• id: participant id
• age is participant age in years
• gender: (m)ale, (f)emale, (d)iverse
• x and y are the sampled coordinates on the grid
• chosen: are the x and y coordinates of the chosen tile
• z is the reward obtained from the chosen tile, normalized to the range 0-1. Re-clicked tiles could show small variations in the observed color (i.e., underlying reward) due to normally distributed noise,\(\epsilon∼N(0,1)\).
• z_scaled is the observed outcome (reward), scaled in each round to a randomly drawn maximum value in the range of 70% to 90% of the highest reward value
• trial is the trial number (0-25), with 0 corresponding to the initially revealed random tile, i.e. trial 1 is the first choice
• round is the round number (1 through 10), with 1=practice round (not analyzed) and 10=bonus round (analyzed only for bonus round judgments)
• distance is the Manhattan distance between consecutive clicks. NA for trial 0, the initially revealed random tile
• type_choice categorizes consecutive clicks as “repeat” (clicking the same tile as in the previous round), “near” (clicking a directly neighboring tile, i.e. distance=1), and “far” (clicking a tile with distance > 1). NA for trial 0, i.e., the initially revealed random tile.
• previous_reward is the reward z obtained on the previous step. NA for trial 0, i.e., the initially revealed random tile.
• last_ldopa: time of the last L-Dopa dose (HH:MM)
• next_ldopa: scheduled time of the next L-Dopa dose (HH:MM)
• time_exp: time of the experiment (HH:MM)
• time_since_ldopa: time since last L-Dopa (in minutes)

File modelFits.csv contains the results of the computational model simulations (GP-UCB model and lesioned variants).

########################################################
# get behavioral data
########################################################
dat_gridsearch <- read_csv("data/data_gridsearch_Parkinson.csv", show_col_types = FALSE) %>% 
  mutate(type_choice  = factor(type_choice, levels = c("Repeat", "Near", "Far"))) 

# normalize reward and previous reward
dat_gridsearch$z = dat_gridsearch$z / 50
dat_gridsearch$previous_reward = dat_gridsearch$previous_reward / 50

########################################################
# get subject data
########################################################
dat_sample <- read_delim("data/data_gridsearch_subjects.csv", escape_double = FALSE, trim_ws = TRUE, show_col_types = FALSE) %>% 
  mutate(gender = as.factor(gender),
        group = fct_recode(group,
                       "Control" = "PNP"
                       # "PD+"     = "PD+",
                       # "PD-"     = "PD-"
                       ),
    group = fct_relevel(group, "PD-", "PD+", "Control")
  ) %>% 
  mutate(last_ldopa = if_else(group != "Control", as_hms(last_ldopa), as_hms(NA)),
         next_ldopa = if_else(group != "Control", as_hms(next_ldopa), as_hms(NA)),
         time_exp = if_else(group != "Control", as_hms(time_exp), as_hms(NA))) %>% 
  mutate(time_since_ldopa = as.numeric(time_exp - last_ldopa, unit = "mins"))


dat <- dat_sample %>% 
  left_join(dat_gridsearch, by = "id") %>% 
  arrange(group)

########################################################
# get modeling data
########################################################
modelFits <- read.csv('modelResults/modelFit.csv') # generated by dataProcessing_gridSearchParkinson.R
# length(unique(modelFits$id))

3 Computational Analyses

Complementing the behavioral analyses, we study exploration and exploitation in PD through the lens of a computational model, the Gaussian Process Upper Confidence Bound (GP-UCB) model. This model integrates similarity-based generalization with two distinct exploration mechanisms: directed exploration, which seeks to reduce uncertainty about rewards, and random exploration, which adds stochastic noise to the search process without being directed towards a particular goal (Wu et al., 2018; Wu et al., 2025). In previous research using the same paradigm, this model has provided the best account of human behavior and enabled the decomposition of exploration into distinct mechanisms (Giron et al., 2023; Meder et al., 2021; Schulz et al., 2019; Wu et al., 2018; Wu et al., 2020).

3.1 Gaussian Process Upper Confidence Bound (GP-UCB) Model

The GP-UCB model comprises three components:

1. a learning model, which uses Bayesian inference to generate predictions about the rewards associated with each option (tile),
2. a sampling strategy, which uses reward expectations and associated uncertainty to evaluate how promising each option is, and
3. a choice rule, which converts options’ values into choice probabilities.

Note

Add details

3.1.1 Learning Model

3.1.2 Sampling Strategy

3.1.3 Choice rule

3.1.4 Model parameters

Associated with each model component is a free parameter that we estimate through out-of-sample cross validation. These parameters provide a window into distinct aspects of learning and exploration:

1. The length-scale parameter \(\lambda\) of the RBF kernel captures how strongly a participant generalizes based on the observed evidence, i.e., the rewards obtained from previous choices.
2. The uncertainty bonus \(\beta\) represents to the level of directed exploration, i.e., how much expected rewards are inflated through an “uncertainty bonus”.
3. The temperature parameter \(\tau\) corresponds to the amount of sampling noise, i.e., extent of random exploration.

4 Model comparison

We tested the GP-UCB model in its ability to model learning and predicting each participants’ search and decision-making behavior. To assess the contribution of each component of the model (generalization, uncertainty-directed exploration, and random exploration) we compare the predictive accuracy of the GP-UCB model to model variants where we lesion away each component.

\(\lambda\) lesion model: This model removes the ability to generalize, meaning that all options are learned independently (via Bayesian mean tracker)

\(\beta\) lesion model: No uncertainty-directed exploration (\(\beta=0\)), i.e., options are valued solely based on reward expectations (mean greedy)

\(\tau\) lesion model: Exchanges the softmax choice rule with an \(\epsilon\)-greedy policy as an alternative random exploration mechanism. With probability \(\epsilon\), a random option is selected (each with probability 1/64); with probability 1 − \(\epsilon\), the option with the highest UCB value is chosen. The parameter \(\epsilon\) is estimated for each participant.

All models were fitted using leave-one-round-out cross-validation based on maximum likelihood estimation. Model fits are evaluated using the sum of negative log-likelihoods across all out-of-sample predictions.

Models’ predictive accuracy was assessed using a pseudo-\(R^2\) measure, based on the sum of negative log-likelihoods across all out-of-sample predictions. The summed log loss is compared to a random model, such that \(R^2=0\) corresponds to chance performance and \(R^2=1\) corresponds to theoretically perfect predictions.

\[ R^2 = 1 - \frac{\log \mathcal{L}(M_k)}{\log\mathcal{L}(M_{rand})}, \]

# get subject information
groupDF <-dat_sample %>%
  select(id, age, gender,group,BDI,MMSE,hoehn_yahr,last_ldopa,next_ldopa,time_exp,time_since_ldopa) %>%
  group_by(id) %>%
  slice_head(n = 1) %>%
  arrange(group)
  
# groupDF <- dat %>% 
#   group_by(id) %>%
#   slice(1) %>%
#   ungroup()

# length(unique(dat$id))
# length(unique(groupDF$id))

modelFits <- merge(modelFits, groupDF[,c('id', 'group')], by = "id") # merge to add group 

# kernels <- c("RBF", "BMT") # RBF = Radial Basis Function kernel, BMT= Bayesian Mean Tracker
# acqFuncs <- c("GM", "UCB", "EG") # UCB = Upper Confidence Bound, GM=greedyMean, EG = epsilonGreedy
modelFits <-  modelFits %>%
  mutate(kernel=factor(kernel, levels=c('RBF', 'BMT'), labels=c('GP', 'BMT'))) %>%
  mutate(acq=factor(acq, levels=c('UCB', 'GM','epsilonGreedy'), labels=c('UCB', 'meanGreedy', 'epsilonGreedy')))

modelFits$ModelName = paste(modelFits$kernel, modelFits$acq, sep="-")

# Only include key comparisons
modelFits <- subset(modelFits, ModelName %in% c("GP-UCB", "BMT-UCB", "GP-meanGreedy", "GP-epsilonGreedy" ))
modelFits$ModelName = factor(modelFits$ModelName, levels = c('GP-UCB', 'BMT-UCB', 'GP-meanGreedy', 'GP-epsilonGreedy'))

#Two line name for models
modelFits$shortname <- factor(modelFits$ModelName, labels = c('GP\nUCB', 'lambda\nlesion', 'beta\nlesion', 'tau\nlesion'))
levels(modelFits$shortname) <- c('GP\nUCB', 'lambda\nlesion', 'beta\nlesion', 'tau\nlesion')

# perform frequentist and Bayesian t-tests and make labels for plotting 
comparisons_df <- modelFits %>%
  group_by(group) %>%
  group_modify(~{
    # pariwise t-tests
    comparisons <- list(
      c("GP\nUCB", "lambda\nlesion"),
      c("GP\nUCB", "beta\nlesion"),
      c("GP\nUCB", "tau\nlesion")
    )
    
    t_res <- t_test(R2 ~ shortname, 
                    data = .x, 
                    paired = TRUE,
                    comparisons = comparisons) %>%
      add_xy_position(x = "shortname") %>%
      mutate(
        p.format = case_when(
          p < 0.001 ~ "p<.001",
          TRUE ~ paste0("p=", signif(p, 2))
        )
      )
    
    # compute Bayes Factors BF10
    t_res$BF <- purrr::pmap_dbl(
      list(t_res$group1, t_res$group2),
      function(g1, g2) {
        x1 <- .x$R2[.x$shortname == g1]
        x2 <- .x$R2[.x$shortname == g2]
        bf <- BayesFactor::ttestBF(x = x1, y = x2, paired = TRUE)
        as.numeric(BayesFactor::extractBF(bf)$bf)
      }
    )
    
    t_res
  }) %>% 
  mutate( # make BF label
    BF.format = case_when(
      BF > 100 ~ "BF>100",
      TRUE ~ paste0("BF=", signif(BF, 2))
    )) %>% 
  mutate(plot_label = paste0(p.format, ", ", BF.format)) # make plot label

# get y.positions from first comaprisons
ref_ypos <- comparisons_df %>%
  filter(group == first(levels(modelFits$group))) %>%
  pull(y.position)

# set manually
ref_ypos <- c(0.68, 0.74, 0.8)

comparisons_df <- comparisons_df %>%
  group_by(group) %>%
  mutate(y.position = ref_ypos) %>%
  ungroup()

p_model_comparison <- ggplot(modelFits, aes(x = shortname, y = R2, fill = group, shape = group, color = group)) +
  facet_wrap(~group, nrow = 1) +
  geom_boxplot(alpha = 0.2, width = 0.4, outlier.shape = NA) +  
  geom_jitter(width = 0.1, size = 1.5, alpha = 0.6) +  
  stat_summary(fun = mean, geom = "point", shape = 23, fill = "white", size = 4) +
  scale_y_continuous(name = expression("Predictive accuracy " ~ R^2),
                     breaks = c(0, 0.5, 1))+  
  scale_x_discrete("",
                   labels = c(
      "lambda\nlesion" = "λ\nlesion",
    "beta\nlesion"   = "β\nlesion",
    "tau\nlesion"    = "τ\nlesion"
    # labels = c(
    # "lambda\nlesion" = expression(atop(lambda, lesion)),
    # "beta\nlesion"   = expression(atop(beta, lesion)),
    # "tau\nlesion"    = expression(atop(tau, lesion))
    )) + 
  scale_fill_manual(values = groupcolors) + 
  scale_color_manual(values = groupcolors) + 
  ggtitle("Model comparison: GP-UCB vs. lesioned models") + 
  theme_classic() +
  theme(strip.background = element_blank(),  
        strip.text = element_text(color = "black", size=18),
        legend.position = "none",
        legend.justification = c(0, 1),
        axis.title = element_text(size = 18),
        axis.text = element_text(size = 18),
        plot.title = element_text(size = 24),
        plot.margin = margin(0, 0, 20, 0) # positive bottom margin, otherwise artfecat when putting later together with cowplot
  )

   
p_model_comparison
ggsave("plots/model_comparison.png", p_model_comparison, dpi=300, width = 10, height = 5)
#ggsave("plots/model_comparison.pdf", p_model_comparison, width = 10, height = 5) # ü

# plot for Computational Psychiatry Conference (CPP; Tübingen, July 2025)
# ggboxplot(modelFits, 
#           x = "shortname", 
#           y = "R2",
#           color = "group", palette = groupcolors, fill = "group", alpha = 0.2,
#           add = "jitter", jitter.size = 1, shape = "group",
#           title = "Model comparison: GP-UCB vs. lesioned models") +
#   facet_wrap(~group, nrow = 1) +
#   ylab(bquote(R^2)) +
#   xlab("") +
#   stat_pvalue_manual(
#     filter(comparisons_df, group == "Control"),
#     label = "plot_label",
#     tip.length = 0.01, bracket.size = 0.3, size = 3
#   ) +
#   stat_pvalue_manual(
#     filter(comparisons_df, group == "PD+"),
#     label = "plot_label",
#     tip.length = 0.01, bracket.size = 0.3, size = 3
#   ) +
#   stat_pvalue_manual(
#     filter(comparisons_df, group == "PD-"),
#     label = "plot_label",
#     tip.length = 0.01, bracket.size = 0.3, size = 3
#   ) +
#   theme_classic()+
#   theme(strip.background = element_blank(),  
#         strip.text = element_text(color = "black", size=12),
#         legend.position = "none",
#         plot.title = element_text(size = 24),
#         axis.text= element_text(colour="black", size = 14),
#         axis.title= element_text(colour="black", size = 14)
#   ) 
# 
# ggsave("plots/model_comparison_CPP.png", p_model_comparison_CPP, dpi=300, width = 9, height = 5)

Figure 1: Predictive accuracy of GP-UCB model and lesioned variants.

4.1 Model comparison: Control

• GP-UCB vs. lambda lesion: \(t(33)=3.0\), \(p=.005\), \(d=0.2\), \(BF=7.5\)
• GP-UCB vs. beta lesion: \(t(33)=3.4\), \(p=.002\), \(d=0.2\), \(BF=19\)
• GP-UCB vs. tau lesion: \(t(33)=7.7\), \(p<.001\), \(d=0.6\), \(BF>100\)

4.2 Model comparison: PD+

• GP-UCB vs. lambda lesion: \(t(32)=3.4\), \(p=.002\), \(d=0.4\), \(BF=20\)
• GP-UCB vs. beta lesion: \(t(32)=3.7\), \(p<.001\), \(d=0.4\), \(BF=40\)
• GP-UCB vs. tau lesion: \(t(32)=8.5\), \(p<.001\), \(d=0.9\), \(BF>100\)

4.3 Model comparison: PD-

• GP-UCB vs. lambda lesion: \(t(30)=3.6\), \(p=.001\), \(d=0.7\), \(BF=27\)
• GP-UCB vs. beta lesion: \(t(30)=5.4\), \(p<.001\), \(d=1.1\), \(BF>100\)
• GP-UCB vs. tau lesion: \(t(30)=4.9\), \(p<.001\), \(d=1.0\), \(BF>100\)

4.4 Model-based classification of participants

# classify participants according to model R^2
df_participant_classification <- modelFits %>%
  group_by(id) %>%
  slice_max(order_by = R2, n = 1) %>%
  select(id, group, ModelName, shortname, R2) %>% 
  ungroup() %>% 
  rename(best_ModelName = ModelName,
         best_shortname = shortname,
         best_R2 = R2)

df_counts <- df_participant_classification %>%
  count(group, best_shortname)

df_percent <- df_counts %>%
  group_by(group) %>%
  mutate(
    total_in_group = sum(n),
    percent = round((n / total_in_group) * 100, 1)
  ) %>%
  ungroup()

# add most predictive model for each subject to df modelFits
modelFits <- modelFits %>% 
  left_join(df_participant_classification, by = c("id", "group"))

We classified participants based on which model achieved the highest cross-validated predictive accuracy (highest \(R^2\); ?@fig-participant_classification). In each patient group, the GP-UCB model was the most predictive model for the majority of participants (Control: 55.9%, PD+: 57.6%, PD-: 58.1%).

In total, out of 98 participants, 56 (57.1%) were best described by the GP-UCB model, 22 (22.4%) by the lambda lesion model, 13 (13.3%) by the beta lesion model, and 7 (7.1%) by the tau lesion model. The results suggest that all three components of the GP-UCB model are relevant for predicting participants’ behavior.

# waffle plot
p_classification_participants <- ggplot(
  data = df_counts, 
  aes(fill=best_shortname, values=n)
) +
  geom_waffle(
    color = "white", 
    size = 1, 
    n_rows = 5
  ) +
  facet_wrap(~group, nrow=1) +
  scale_x_discrete(
    expand = c(0,0,0,0)
  ) +
  scale_y_discrete(
    expand = c(0,0,0,0)
  ) +
  ggthemes::scale_fill_tableau(name=NULL) +
  coord_equal() +
  ggtitle ("Model comparison: Participant classification") +
theme_classic() +
  theme(
    legend.title = element_blank(),
    plot.title = element_text(size = 24),
    legend.position = 'right',
    strip.text = element_text(color = "black", size=18),
    legend.text =  element_text(colour="black", size=18),
    text = element_text(colour = "black"),
    strip.background =element_blank(),
    axis.text= element_text(colour="black", size = 18),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.spacing = unit(3, "lines"),
    legend.key.spacing.y = unit(0.4, "cm"),
    plot.margin = margin(-70, 0, -30, 0) # negativ bottom margin, otherwise artfecat when putting later together with cowplot
    )

ggsave("plots/participant_classification.png", p_classification_participants, width = 12, height = 5, dpi=300)

5 Analysis of parameter estimates

df_gpucb_params <- modelFits %>% filter(kernel=='GP' & acq == 'UCB') %>% 
  pivot_longer(c('lambda', 'beta', 'tau'), names_to = 'param', values_to = 'estimate') %>% 
  mutate(param = factor(param, levels = c('lambda', 'beta', 'tau'))) %>% 
  mutate(estimate_log10 = log10(estimate))


df_gpucb_params %>%
  group_by(group, param) %>%
  summarise(
    mean = mean(estimate, na.rm = TRUE),
    median = median(estimate, na.rm = TRUE),
    se = sd(estimate, na.rm = TRUE) / sqrt(n()),
    ci_lower = mean - 1.96 * se,
    ci_upper = mean + 1.96 * se,
    .groups = "drop"
  ) %>%
  mutate(
    summary = sprintf("M = %.2f [%.2f, %.2f], Mdn = %.2f", mean, ci_lower, ci_upper, median)
  ) %>%
  select(group, param, summary) %>%
  pivot_wider(names_from = param, values_from = summary) %>%
  kable(caption = "Mean (95% CI) and median parameter estimates by group", format = "html") %>% 
  kable_styling("striped", full_width = FALSE)

Table 1: Mean (95% CI) and median parameter estimates of the GP-UCB model by group.

Mean (95% CI) and median parameter estimates by group

group	lambda	beta	tau
PD-	M = 0.54 [0.46, 0.63], Mdn = 0.46	M = 10.22 [3.77, 16.68], Mdn = 0.55	M = 1.35 [0.20, 2.50], Mdn = 0.05
PD+	M = 0.56 [0.50, 0.62], Mdn = 0.53	M = 1.99 [-0.22, 4.20], Mdn = 0.37	M = 0.42 [-0.11, 0.94], Mdn = 0.04
Control	M = 0.65 [0.57, 0.73], Mdn = 0.59	M = 2.36 [-0.64, 5.35], Mdn = 0.37	M = 0.65 [-0.40, 1.71], Mdn = 0.04

# separate plots
# generalization lambda
p_gpucb_params_lambda <- 
  ggplot(subset(df_gpucb_params, param == "lambda"), aes(x = group, y = estimate, fill = group, shape = group, color = group)) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "darkred") + # true lambda
  scale_y_continuous(name = "Estimate", breaks = c(0, 0.5, 1, 2), labels = c("0", "0.5", "1", "2"), limits = c(0, 1.45)) +
  geom_boxplot(alpha = 0.2, width = 0.4, outlier.shape = NA) +  
  geom_jitter(width = 0.1, size = 1.5, alpha = 0.6) +  
  stat_summary(fun = mean, geom = "point", shape = 23, fill = "white", size = 4) +
  scale_color_manual(values = groupcolors) +
  scale_fill_manual(values = groupcolors) +
  scale_x_discrete("") + 
  # ggtitle("GP-UCB parameter estimates: Group differences") + 
  ggtitle(expression("Generalization " * lambda)) + 
  theme_classic() +
  theme(legend.position = "none",
        legend.justification = c(0, 1),
        axis.title = element_text(size = 14),
        axis.text = element_text(size = 14),
        plot.title = element_text(size = 18, hjust = 0.5),
        plot.margin = margin(10, 0, 10, 0) 
  )

# p_gpucb_params_lambda  
# ggsave("plots/GP-p_gpucb_params_lambda.png", p_gpucb_params_lambda, dpi=300, width = 4, height = 4)

# exploration bonus beta
p_gpucb_params_beta <- 
  ggplot(subset(df_gpucb_params, param == "beta"), aes(x = group, y = estimate, fill = group, shape = group, color = group)) +
  scale_y_log10(name = " ", breaks = c(0.01, 0.1, 1, 10), labels = c("0.01", "0.1", "1", "10"), limits = c(0.005, 55)) +
  geom_boxplot(alpha = 0.2, width = 0.4, outlier.shape = NA) +  
  geom_jitter(width = 0.1, size = 1.5, alpha = 0.6) +  
  stat_summary(fun = mean, geom = "point", shape = 23, fill = "white", size = 4) +
  scale_color_manual(values = groupcolors) +
  scale_fill_manual(values = groupcolors) +
  scale_x_discrete("") + 
  # ggtitle("GP-UCB parameter estimates: Group differences") + 
  ggtitle(expression("Exploration bonus " * beta)) + 
  theme_classic() +
  theme(legend.position = "none",
        legend.justification = c(0, 1),
        axis.title = element_text(size = 14),
        axis.text = element_text(size = 14),
        plot.title = element_text(size = 18, hjust = 0.5),
        plot.margin = margin(10, 0, 10, 0) 
  )

# p_gpucb_params_beta  
# ggsave("plots/GP-p_gpucb_params_beta.png", p_gpucb_params_beta, dpi=300, width = 4, height = 4)

# random exploration temperature tau 
p_gpucb_params_tau <- 
  ggplot(subset(df_gpucb_params, param == "tau"), aes(x = group, y = estimate, fill = group, shape = group, color = group)) +
  scale_y_log10(name = " ", breaks = c(0.01, 0.1, 1, 10), labels = c("0.01", "0.1", "1", "10"), limits = c(0.005, 55)) +
  geom_boxplot(alpha = 0.2, width = 0.4, outlier.shape = NA) +  
  geom_jitter(width = 0.1, size = 1.5, alpha = 0.6) +  
  stat_summary(fun = mean, geom = "point", shape = 23, fill = "white", size = 4) +
  scale_color_manual(values = groupcolors) +
  scale_fill_manual(values = groupcolors) +
  scale_x_discrete("") + 
  # ggtitle("GP-UCB parameter estimates: Group differences") + 
  ggtitle(expression("Random exploration " * tau)) + 
  theme_classic() +
  theme(legend.position = "none",
        legend.justification = c(0, 1),
        axis.title = element_text(size = 14),
        axis.text = element_text(size = 14),
        plot.title = element_text(size = 18, hjust = 0.5),
        plot.margin = margin(10, 0, 10, 0) 
  )

# p_gpucb_params_tau  
# ggsave("plots/GP-p_gpucb_params_tau.png", p_gpucb_params_tau, dpi=300, width = 4, height = 4)

# out together and add title
p_gpucb_params <- grid.arrange(
  grobs = list(p_gpucb_params_lambda, p_gpucb_params_beta, p_gpucb_params_tau),
  nrow = 1,
  top = textGrob(
    "         GP-UCB parameters: Group differences", 
    x = 0,            
    hjust = 0,        
    gp = gpar(fontsize = 24)
  ),
  padding = unit(0.5, "lines") 
)

ggsave("plots/gpucb_params.png", p_gpucb_params, dpi = 300, width = 12, height = 4)

Figure 2: Parameter estimates of GP-UCB model, estimated through leave-one-round-out cross validation. Each dot is one participant.

# plot for Computational Psychiatry Conference (CPP; Tübingen, July 2025)
# 
# # Define your comparisons
# comparisons <- list(c("PD+", "PD-"), c("Control", "PD+"), c("Control", "PD-"))
# 
# # Extract function for p and BF from ttestPretty output
# # TO DO: Cumbersome via 
# extract_p_and_bf <- function(tt_string) {
#   matches <- stringr::str_match_all(tt_string, "\\$p=([^$]+)\\$|\\$BF=([^$]+)\\$")
#   flat <- unlist(matches)
#   
#   raw_vals <- flat[!is.na(flat) & grepl("^\\.?\\d+", flat)]
#   
#   # Convert to numeric and round to 2 decimal places
#   nums <- signif(as.numeric(raw_vals), 2)
#   
#   # Format with 2 decimal digits (or scientific if very small/large)
#   p_fmt <- formatC(nums[1], digits = 2, format = "f")
#   bf_fmt <- formatC(nums[2], digits = 2, format = "f")
#   
#   paste0("p=", p_fmt, ", BF=", bf_fmt)
# }
# 
# # Loop over each param and each comparison
# comparisons_df <- df_gpucb_params %>%
#   group_by(param) %>%
#   group_modify(~{
#     comparisons <- list(
#       c("Control", "PD+"),
#       c("Control", "PD-"),
#       c("PD+", "PD-")
#     )
#     
#     # For each pairwise group comparison
#     res <- purrr::map_dfr(comparisons, function(groups) {
#       g1 <- groups[1]
#       g2 <- groups[2]
#       
#       x1 <- .x$estimate[.x$group == g1]
#       x2 <- .x$estimate[.x$group == g2]
#       
#       if (length(x1) < 2 || length(x2) < 2) {
#         return(tibble(
#           group1 = g1,
#           group2 = g2,
#           p = NA_real_,
#           BF = NA_real_,
#           y.position = NA_real_
#         ))
#       }
#       
#       # Frequentist test
#       t_res <- t.test(x1, x2, paired = FALSE, var.equal = TRUE)
#       p_val <- t_res$p.value
#       
#       # Bayes Factor
#       bf <- BayesFactor::ttestBF(x = x1, y = x2, paired = FALSE)
#       bf_val <- as.numeric(BayesFactor::extractBF(bf)$bf)
#       
#       # y-position (max value in current param group × offset)
#       y_max <- max(.x$estimate, na.rm = TRUE)
#       y_pos <- y_max * runif(1, 1.05, 1.15)
#       
#       tibble(
#         group1 = g1,
#         group2 = g2,
#         p = p_val,
#         BF = bf_val,
#         y.position = y_pos
#       )
#     })
#     
#     res
#   }) %>%
#   ungroup() %>%
#   mutate(
#     p.format = case_when(
#       p < 0.001 ~ "p<.001",
#       is.na(p) ~ NA_character_,
#       TRUE ~ paste0("p=", signif(p, 2))
#     ),
#     BF.format = case_when(
#       is.na(BF) ~ NA_character_,
#       BF > 100 ~ "BF>100",
#       TRUE ~ paste0("BF=", signif(BF, 2))
#     ),
#     plot_label = paste0(p.format, ", ", BF.format)
#   )
# 
# comparisons_df$y.position <- rep(c(1.8, 2.6, 2.2), 3)
# 
# p_GP_UCB_params_CPP <-  
#   ggboxplot(df_gpucb_params, 
#           x = "group", 
#           y = "estimate",
#           color = "group", palette =groupcolors, fill = "group", alpha = 0.2,
#           add = "jitter", jitter.size = 0.5, shape = "group", title = "GP-UCB parameter estimates") +
#   facet_wrap(~param, nrow = 1) +
#   scale_y_log10(breaks = c(0.01, 0.1, 1, 10, 100), labels = c("0.01", "0.1", "1", "10", "100"), expand = expansion(mult = c(0.1, 0.15))  ) +
#   # scale_y_log10( expand = expansion(mult = c(0, 0.1))  ) +
#   # coord_cartesian(ylim = c(0.01,260)) +
#   ylab("Estimate (log scale)") +
#   xlab("") +
#     # ignore p values because of log; only done to position brackets correctly
#  # stat_compare_means(comparisons = list( c("Control", "PD+"), c("PD+", "PD-"), c("Control", "PD-")  ),
#  #                     paired = F,
#  #                     method = "t.test",
#  #                     # label = "p.format",
#  #                     aes(label = paste0("p = ", after_stat(p.format)))
#  #                     #aes(label = paste0(" "))
#  #  ) +
#     stat_pvalue_manual(
#     filter(comparisons_df, param == "lambda"),
#     label = "plot_label",
#     tip.length = 0.01, bracket.size = 0.3, size = 3
#   ) +
#     stat_pvalue_manual(
#     filter(comparisons_df, param == "beta"),
#     label = "plot_label",
#     tip.length = 0.01, bracket.size = 0.3, size = 3
#   ) +
#     stat_pvalue_manual(
#     filter(comparisons_df, param == "tau"),
#     label = "plot_label",
#     tip.length = 0.01, bracket.size = 0.3, size = 3
#   ) +
#   stat_summary(fun = mean, geom="point", shape = 23, fill = "white", size=2) +
#   theme_classic() +
#   theme(strip.background = element_blank(),  
#         strip.text = element_text(color = "black", size=14),
#         legend.position = "none",
#         plot.title = element_text(size = 20),
#         axis.text= element_text(colour="black", size = 12),
#         axis.title= element_text(colour="black", size = 12),
#         panel.spacing = unit(3, "lines") 
#   )
# 
# ggsave("plots/GP-UCB_params_CPP.png", p_GP_UCB_params_CPP, dpi = 300, width = 8, height = 5)
# ggsave("plots/GP-UCB_params_CPP.pdf", p_GP_UCB_params_CPP, width = 8, height = 5)

To better understand the mechanisms underlying the observed behavioral differences, we analyzed the parameters of the Gaussian Process Upper Confidence Bound (GP-UCB) model (Figure 2).

5.0.1 Generalization \(\lambda\)

The parameter \(\lambda\) represents the length-scale in the RBF kernel, which governs the amount of generalization, i.e., to what extent participants assume a spatial correlation between options (higher \(\lambda\) = stronger generalization). Overall, the amount of generalization was very similar between groups.

• Control vs. PD+: \(U=678\), \(p=.145\), \(r_{ au}=.15\), \(BF=.65\)
• Control vs. PD-: \(U=731\), \(p=.007\), \(r_{ au}=.28\), \(BF=3.9\)
• PD+ vs. PD-: \(U=626\), \(p=.126\), \(r_{ au}=.16\), \(BF=.44\)

5.0.2 Exploration bonus \(\beta\)

The parameter \(\beta\) represents the uncertainty bonus, i.e. how much expected rewards are positively inflated by their uncertainty (higher \(\beta\) = more uncertainty-directed exploration). Controls and PD+ patients on medication did not differ, and both groups had lower beta estimates than the dopamine-depleted patients in the PD− group. These differences suggest that levodopa medication modulated the amount of uncertainty-directed exploration by restoring beta to levels comparable to those observed in controls without PD. This aligns with findings from a restless bandit paradigm, where L-Dopa reduced the amount of directed exploration in healthy volunteers, while the level of random exploration remained unaffected (Chakroun et al., 2020).

• Control vs. PD+: \(U=480\), \(p=.315\), \(r_{ au}=-.10\), \(BF=.48\)
• Control vs. PD-: \(U=188\), \(p<.001\), \(r_{ au}=-.46\), \(BF=81\)
• PD+ vs. PD-: \(U=220\), \(p<.001\), \(r_{ au}=-.41\), \(BF=25\)

5.0.3 Random exploration \(\tau\)

The parameter \(\tau\) represents the amount of decision noise, i.e. stochastic variability in the softmax decision rule (lower \(\tau\) = more decision noise, i.e. more uniform distribution; conversely, \(\tau \rightarrow \infty \quad \Rightarrow \quad \text{argmax (greedy)}\)). There were no group differences in rge temperature paramter \(\tau\), indicating comparable amounts of random exploration regardless of group.

• Control vs. PD+: \(U=572\), \(p=.896\), \(r_{ au}=.01\), \(BF=.25\)
• Control vs. PD-: \(U=500\), \(p=.730\), \(r_{ au}=-.04\), \(BF=.27\)
• PD+ vs. PD-: \(U=470\), \(p=.584\), \(r_{ au}=-.06\), \(BF=.28\)

6 Relations of model parameters to performance

We assessed the correlation (Kendall’s tau, because it’s invariant against log transformation) of GP-UCB parameter estimates with performance (mean reward).

# mean reward per subject across all trials and rounds (practice and bonus round excluded)
df_mean_reward_subject <- dat %>% 
  filter(trial != 0 & round %in% 2:9) %>% # exclude first (randomly revealed) tile and practice round and bonus round
  group_by(id) %>% 
  summarise(group = first(group),
            sum_reward = sum(z),
            mean_reward = mean(z), 
            sd_reward = sd(z)) 

df_params_performance <- df_gpucb_params %>% 
  left_join(df_mean_reward_subject, by = c("id", "group"))

df_params_performance_wide <- df_gpucb_params %>% 
  pivot_wider(names_from = param, values_from = estimate ) %>% 
  left_join(df_mean_reward_subject, by = c("id", "group"))

The amount of generalization was positively related with obtained rewards, showing that participants who successfully learned about the spatially correlation of rewards performed better. The uncertainty bonus \(\beta\) was negatively correlated with performance, demonstrating that an overreliance on exploration impairs efficient reward accumulation. The amount of random temperature \(\tau\) was not related to obtained rewards.

# plot correlation lambda and reward 
p_lambda_reward <- 
  ggplot(subset(df_params_performance, param == "lambda"), aes(x = estimate, y = mean_reward)) +
  geom_point(aes(color = group, shape = group, fill = group)) +
  # geom_smooth(method = "lm", color = "black", linetype = "dashed", fill = "lightgray", se = TRUE) + # one regression line
  geom_smooth(aes(color = group, fill = group), method = "lm", linetype = "dashed", alpha = 0.2, se = TRUE) +  # one regrssion line per group
  stat_cor(method = "kendall", cor.coef.name = expression(r[tau]), label.x.npc = "left",  label.y = 0.4,  size=5, p.accuracy=0.001) +
  labs(
    title = expression("Generalization " * lambda),
    x = "Estimate (log)",
    y = "Mean normalized reward"
  ) +
  # scale_x_log10(name = "Estimate (log scale)", breaks = c(0.1, 1, 10), labels = c("0.1", "1", "10"), limits = c(0.1, 2)) +
  scale_x_continuous(name = "Estimate", breaks = c(0, 1, 2), labels = c("0", "1", "2"), limits = c(0, 1.45)) +
  scale_y_continuous(breaks = seq(0, 1, .1), limits = c(0.4, 0.85)) +
  scale_fill_manual(values = groupcolors) + 
  scale_color_manual(values = groupcolors) +
  theme_classic() +
  theme(
    plot.title = element_text(hjust = 0.5, size=18),
    legend.title = element_blank(),
    axis.title = element_text(color = "black", size = 14),
    axis.text = element_text(color = "black", size = 14),
    legend.position = c(0.12, 0.9), #bottom: 0.15 top:0.9
    legend.text = element_text(size = 14),
    legend.key.size = unit(0.8, "lines"),
    legend.margin = margin(0, 0, 0, 0),       
  legend.box.margin = margin(0, 0, 0, 0),
  plot.margin = margin(10, 0, 0, 0), 
  legend.key = element_rect(fill = NA, colour = NA)

  )

# ggsave("plots/cor_lambda_reward.png", p_lambda_reward, dpi = 300, width = 8, height = 5)

# plot correlation between exploration bonus beta (log10) and reward  
p_beta_reward <-
  ggplot(subset(df_params_performance, param == "beta"), aes(x = estimate, y = mean_reward)) +
   geom_point(aes(color = group, shape = group, fill = group)) +
  # geom_smooth(method = "lm", color = "black", linetype = "dashed", fill = "lightgray", se = TRUE) + # one regression line
  geom_smooth(aes(color = group, fill = group), method = "lm", linetype = "dashed", alpha = 0.2, se = TRUE) +  # one regrssion line per group
  stat_cor(method = "kendall", cor.coef.name = expression(r[tau]), label.x.npc = "left",  label.y = 0.4,  p.accuracy=0.001, size=5) + #label.x = 0.00,  label.y = 0.4,
  labs(
    title = expression("Exploration bonus " * beta),
    x = "Estimate (log)",
    y = " "
  ) +
  scale_x_log10(name = "Estimate (log scale)",  breaks = c(0.01, 0.1, 1, 10), labels = c("0.01","0.1", "1", "10"), limits = c(0.005, 70)) +
  scale_y_continuous(breaks = seq(0, 1, .1), limits = c(0.4, 0.85)) +
  scale_fill_manual(values = groupcolors) + 
  scale_color_manual(values = groupcolors) +
  theme_classic() +
  theme(
    plot.title = element_text(hjust = 0.5, size=18),
    legend.title = element_blank(),
    axis.title = element_text(color = "black", size = 14),
    axis.text = element_text(color = "black", size = 14),
    legend.position = "none",
    plot.margin = margin(10, 0, 0, 0)
  )

# ggsave("plots/cor_beta_reward.png", p_beta_reward, dpi = 300, width = 8, height = 5)


# plot correlation between amount of random exploration (temperature tau of softmax choice rule) and reward  
p_tau_reward <-
  ggplot(subset(df_params_performance, param == "tau"), aes(x = estimate, y = mean_reward)) +
  geom_point(aes(color = group, shape = group, fill = group)) +
  # geom_smooth(method = "lm", color = "black", linetype = "dashed", fill = "lightgray", se = TRUE) + # one regression line
  geom_smooth(aes(color = group, fill = group), method = "lm", linetype = "dashed", alpha = 0.2, se = TRUE) +  # one regrssion line per group
  stat_cor(method = "kendall", cor.coef.name = expression(r[tau]), label.x.npc = "left",  label.y = 0.4, p.accuracy=0.01, size=5) +
  labs(
    title = expression("Random exploration " * tau),
    x = "Estimate (log)",
    y = " "
  ) +
  scale_x_log10(name = "Estimate (log scale)",  breaks = c(0.01, 0.1, 1, 10), labels = c("0.01","0.1", "1", "10"), limits = c(0.005, 20)) +
  scale_y_continuous(breaks = seq(0, 1, .1), limits = c(0.4, 0.85)) +
  scale_fill_manual(values = groupcolors) + 
  scale_color_manual(values = groupcolors) +
  theme_classic() +
  theme(
    plot.title = element_text(hjust = 0.5, size=18),
    legend.title = element_blank(),
    axis.title = element_text(color = "black", size = 14),
    axis.text = element_text(color = "black", size = 14),
    legend.position = "none",
    plot.margin = margin(10, 0, 0, 0)
  )

# ggsave("plots/cor_tau_reward.png", p_tau_reward, dpi = 300, width = 8, height = 5)

# put plots together and add title
p_parameters_reward <- grid.arrange(
  grobs = list(p_lambda_reward, p_beta_reward, p_tau_reward),
  nrow = 1,
  top = textGrob(
    "         GP-UCB parameters: Relations to performance", 
    x = 0,            
    hjust = 0,        
    gp = gpar(fontsize = 24)
  ),
  padding = unit(0.5, "lines") 
)

p_parameters_reward

TableGrob (2 x 3) "arrange": 4 grobs
  z     cells    name                 grob
1 1 (2-2,1-1) arrange       gtable[layout]
2 2 (2-2,2-2) arrange       gtable[layout]
3 3 (2-2,3-3) arrange       gtable[layout]
4 4 (1-1,1-3) arrange text[GRID.text.1160]

ggsave("plots/parameters_reward.png", p_parameters_reward, dpi = 300, width = 13, height = 4)
 
# # plot correlation between parameter estimates and mean reward, with inset for smaller estimate range 
# main_plot <- ggscatter(df_params_performance, x = "estimate", y = "mean_reward",
#                        add = "reg.line",  
#                        add.params = list(color = "darkred", fill = "lightgray"), 
#                        conf.int = TRUE 
# ) +
#   facet_wrap(~param, scales = "free_x") +
#   stat_cor(method = "pearson", label.x = c(0,3,3), label.y = 45) +
#   ggtitle("Correlation between GP-UCB parameters and reward") +
#   scale_y_continuous("Mean reward", breaks = seq(0,45,10)) +
#   xlab("Estimate") +
#   theme_classic() +
#   theme(strip.background = element_blank(),  
#         strip.text = element_text(color = "black", size=12),
#         legend.title = element_blank(),
#         axis.title = element_text(color = "black", size=14),
#         axis.text = element_text(color = "black", size=14))
# 
# # Inset-Plot für beta (nur Werte 0-1)
# inset_beta <- 
#   ggscatter(df_params_performance %>% filter(param == "beta" & estimate > 0 & estimate <= 1), 
#             x = "estimate", y = "mean_reward",
#             add = "reg.line",  
#             add.params = list(color = "darkred", fill = "lightgray"), 
#             conf.int = TRUE 
#   ) +
#   stat_cor(method = "pearson", label.x = 0.1, label.y = 42, size = 3) +
#   scale_x_continuous( breaks = c(0,0.25, 0.5, 0.75), labels = c("0.0", "0.25", "0.5", "0.75")) +# breaks = seq(0,1,0.1),
#   theme_classic() +
#   theme(axis.title = element_blank(),
#         strip.background = element_blank(), 
#         strip.text = element_blank(), 
#         legend.position = "none")
# 
# # Inset-Plot für tau (nur Werte 0-1)
# inset_tau <- 
#   ggscatter(df_params_performance %>% filter(param == "tau" & estimate > 0 & estimate <= 1), 
#             x = "estimate", y = "mean_reward",
#             add = "reg.line",  
#             add.params = list(color = "darkred", fill = "lightgray"), 
#             conf.int = TRUE 
#   ) +
#   stat_cor(method = "pearson", label.x = 0.05, label.y = 42, size = 3) +
#   scale_x_continuous(breaks = seq(0,1,0.1)) +
#   theme_classic() +
#   theme(axis.title = element_blank(),strip.background = element_blank(), strip.text = element_blank(), legend.position = "none")
# 
# 
# # Hauptplot mit Insets für beta und tau
# p_GP_UCB_params_cor_reward_inset <- 
#   ggdraw(main_plot) +
#   draw_plot(inset_beta, x = 0.5, y = 0.45, width = 0.15, height = 0.35) +
#   draw_plot(inset_tau, x = 0.8, y = 0.45, width = 0.15, height = 0.35)
# 
# p_GP_UCB_params_cor_reward_inset
# 
# ggsave("plots/GP-UCB_params_cor_reward.png", p_GP_UCB_params_cor_reward_inset, width = 12, height= 4, dpi=300)

Figure 3: Correlation of GP-UCB parameters with obtained mean reward across all trials and rounds. Each dot is one participant. The insets show the correlations for a restricted parameter range from 0 to 1.

6.1 Generalization \(\lambda\)

Overall, the extent of generalization was positively related to performance, suggesting that participants who stronger generalized obtained more rewards:

• Overall: \(r_{ au}=.26\), \(p<.001\), \(BF>100\)

Analysis of parameter estimates on the group level showed that this overall relation was primarily driven by PD+ patients, who showed a strong relation, whereas there was no relation in controls or PD- patients:

• Control: \(r_{ au}=.13\), \(p=.288\), \(BF=.39\)
• PD+: \(r_{ au}=.45\), \(p<.001\), \(BF>100\)
• PD-: \(r_{ au}=-.01\), \(p=.973\), \(BF=.23\)

6.2 Exploration bonus \(\beta\)

The exploration bonus \(\beta\) driving uncertainty-directed correlation was negatively related to performance, suggesting that participants who explore too much at the cost of exploiting known high-value options achieve lower performance:

• Overall: \(r_{ au}=-.59\), \(p<.001\), \(BF>100\)

Analysis of parameter estimates on the group level showed that this overall relation was primarily driven by PD+ patients, who showed a strong relation, whereas there was no relation in controls or PD- patients:

• Control: \(r_{ au}=-.43\), \(p<.001\), \(BF>100\)
• PD+: \(r_{ au}=-.61\), \(p<.001\), \(BF>100\)
• PD-: \(r_{ au}=-.60\), \(p<.001\), \(BF>100\)

6.3 Random exploration \(\tau\)

The temperature parameter of the softmax choice rule \(\tau\), representig random exploration, was not related to performance, suggesting that participants who explore too much at the cost of exploiting known high-value options achieve lower performance:

• Overall: \(r_{ au}=-.07\), \(p=.308\), \(BF=.22\)

Analysis of parameter estimates on the group level showed that this overall relation was primarily driven by PD+ patients, who showed a strong relation, whereas there was no relation in controls or PD- patients:

• Control: \(r_{ au}=-.02\), \(p=.860\), \(BF=.23\)
• PD+: \(r_{ au}=.09\), \(p=.451\), \(BF=.30\)
• PD-: \(r_{ au}=-.23\), \(p=.077\), \(BF=1.1\)

6.4 Regression: Perfomance and model parameters

We also performed a regression analysis where we included all model paramaters together with group as predictors.

df_params_performance_wide <- 
  df_params_performance %>% 
    filter(ModelName == "GP-UCB") %>% 
    select(id, group, mean_reward, param, estimate_log10) %>%
    distinct(id, group, param, .keep_all = TRUE) %>%   # drop duplicates if any
  pivot_wider(names_from = param, values_from = estimate_log10) %>%
  drop_na(beta, tau, lambda)   

lm_performance_parameters_log <- lm(mean_reward ~ group * (lambda + beta + tau), data =  df_params_performance_wide)

tab_model(lm_performance_parameters_log)
# res.table <- as.data.frame(coef(summary(lm_performance_parameters_log)))

# check models
# library(performance)
# check_model(lm_performance_parameters_log)
# 
# df_params_performance_wide2 <- 
#   df_params_performance %>% 
#     filter(ModelName == "GP-UCB") %>% 
#     select(id, group, mean_reward, param, estimate) %>%
#     distinct(id, group, param, .keep_all = TRUE) %>%   # drop duplicates if any
#   pivot_wider(names_from = param, values_from = estimate) %>%
#   drop_na(beta, tau, lambda)   
# 
# 
# lm_performance_parameters <- lm(mean_reward ~ group * (lambda + beta + tau), data =  df_params_performance_wide2)
# 
# tab_model(lm_performance_parameters, title = "Regression results: Performance (obtained rewards) as function of group and model parameters.")
# 
# res.table <- as.data.frame(coef(summary(lm_performance_parameters)))
# 
# check_model(lm_performance_parameters)

Table 2: Regression results: Performance (obtained rewards) as function of group and model parameters (log scale).

	mean reward
Predictors	Estimates	CI	p
(Intercept)	0.54	0.50 – 0.58	<0.001
group [PD+]	0.10	0.03 – 0.17	0.009
group [Control]	0.08	0.03 – 0.14	0.005
lambda	0.02	-0.20 – 0.24	0.858
beta	-0.02	-0.06 – 0.01	0.206
tau	-0.01	-0.05 – 0.03	0.608
group [PD+] × lambda	0.07	-0.26 – 0.41	0.675
group [Control] × lambda	-0.10	-0.51 – 0.31	0.621
group [PD+] × beta	-0.05	-0.14 – 0.05	0.345
group [Control] × beta	-0.07	-0.15 – 0.01	0.078
group [PD+] × tau	0.03	-0.04 – 0.09	0.447
group [Control] × tau	0.03	-0.06 – 0.12	0.519
Observations	98
R2 / R2 adjusted	0.635 / 0.588

6.5 Model simulations

To evaluate how well different parameter settings balance exploration and exploitation, we conducted simulations with the GP-UCB model. In these simulations, we fixed the value of \(\lambda\) at 1, corresponding to the true amount of correlation in the used environments, and systematically varied the amount of random exploration (\(\tau\)) and the size of the uncertainty bonus (\(\beta\)). For each parameter we defined used equally log-spaced values, and then simulated 100 learners searching for rewards. Environments were sampled (with replacement) from the set of 40 environments used in the empirical study.

# simulation results
sim = read.csv('modelResults/simulatedModels_local_lambda_1.csv')
#sim = read.csv('modelResults/simulatedModels_local_lambda_0_5.csv')

# normalize reward
sim$meanReward = sim$mu / 50

sim_means <-  sim %>% 
  group_by(tau,beta) %>% 
  summarise(meanReward = mean(mu)) %>% 
  mutate(meanReward = meanReward/50) %>% 
  mutate(beta_log10 = log10(beta),
         tau_log10 = log10(tau)
  )


# get mean parameter estimates by group
marker <- df_gpucb_params %>%
  group_by(group, param) %>%
  summarise(#mean = mean(estimate, na.rm = TRUE), 
            median = median(estimate, na.rm = TRUE),
                             .groups = "drop") %>%
  filter(param %in% c("beta", "tau")) %>%
  pivot_wider(names_from = param, values_from = median) %>%
  mutate(beta_log10 = log10(beta), tau_log10 = log10(tau))

# get median parameter estimates by group
marker2 <-
  df_gpucb_params %>%
  select(id, group, param, estimate) %>%
  filter(param %in% c("beta", "tau")) %>%
  pivot_wider(names_from = param, values_from = estimate) %>%
  mutate(beta_log10 = log10(beta), tau_log10 = log10(tau))

# tick positions in the original scale
bx <- c(0.001, 0.01, 0.1, 1, 10, 50)
by <- c(0.001, 0.01, 0.1, 1, 10, 20)

# Control group
p_model_simulation_params_control <- 
ggplot(sim_means, aes(x = beta_log10, y = tau_log10, fill = meanReward)) +
  geom_raster() +
  scale_x_continuous(breaks = log10(bx), labels = bx, expand = c(0,0)) +
  scale_y_continuous(breaks = log10(by), labels = by, expand = c(0,0)) +
  labs(x = expression(paste('Exploration bonus ', beta)),
       y = expression(paste('Random exploration ', tau))) +
  # scale_fill_gradientn(colours = colorspace::sequential_hcl(500, "plasma"),name = "Normalized\nreward") +
  # scale_fill_viridis_c(option = "plasma", name = "Normalized\nreward") +
  scale_fill_gradientn(colours = cet_pal(5, name = "l7"), name = "Normalized\nreward") + #l3 l6 l16 i5 "rainbow"
  #coord_cartesian(xlim = c(-3,0.5), ylim=c(-3, 0.5)) +
  geom_jitter( # add individual points
    data = subset(marker2, group == "Control"),
    aes(x = beta_log10, y = tau_log10, shape = group, colour = group),
    inherit.aes = FALSE,  
    size = 3,
    fill = "#7570b3",
    shape = 21,
    stroke = 0.3,
    color = "white",
    width = 0.1
  ) +
 geom_point( # add median
    data = subset(marker, group == "Control"),
    aes(x = beta_log10, y = tau_log10, colour = group),
    inherit.aes = FALSE,
    shape = 21,        
    size = 5,          
    stroke = 1.5,
    fill = "#7570b3",
    color = "white"
  ) +
  ggtitle("Control") +
  theme_classic() +
  theme(strip.background = element_blank(),  
        strip.text = element_text(color = "black", size=18),
        legend.position = "none",
        #        legend.justification = c(0, 1),
         
        axis.title = element_text(size = 18),
        axis.text = element_text(size = 18),
        plot.title = element_text(size = 18, hjust = 0.5),
        axis.title.y = element_blank()
  )

#p_model_simulation_params_control

# PD- group
p_model_simulation_params_pd_off <- 
ggplot(sim_means, aes(x = beta_log10, y = tau_log10, fill = meanReward)) +
  geom_raster() +
  scale_x_continuous(breaks = log10(bx), labels = bx, expand = c(0,0)) +
  scale_y_continuous(breaks = log10(by), labels = by, expand = c(0,0)) +
  labs(x = expression(paste('Exploration bonus ', beta)),
       y = expression(paste('Random exploration ', tau))) +
  # scale_fill_gradientn(colours = colorspace::sequential_hcl(500, "plasma"),name = "Normalized\nreward") +
  # scale_fill_viridis_c(option = "plasma", name = "Normalized\nreward") +
   scale_fill_gradientn(colours = cet_pal(5, name = "l7"), name = "Normalized\nreward") + #l3 l16 i5 "rainbow"
  coord_cartesian(xlim = c(-3,0.5), ylim=c(-3, 0.5)) +
  geom_jitter( # add individual points
    data = subset(marker2, group == "PD-"),
    aes(x = beta_log10, y = tau_log10, shape = group, colour = group),
    inherit.aes = FALSE,  
    size = 3,
    color = "white",
    fill = "#d95f02",
    shape = 22,
    stroke = 0.3,
    width = 0.1
  ) +
 geom_point( # add median
    data = subset(marker, group == "PD-"),
    aes(x = beta_log10, y = tau_log10, colour = group),
    inherit.aes = FALSE,
    shape = 22,        
    size = 5,          
    stroke = 1.5,
    fill = "#d95f02",
    color = "white"
  ) +
  ggtitle("PD-") +
  theme_classic() +
  theme(strip.background = element_blank(),  
        strip.text = element_text(color = "black", size=18),
        legend.position = "none",
        #        legend.justification = c(0, 1),
        axis.title = element_text(size = 18),
        axis.text = element_text(size = 18),
        plot.title = element_text(size = 18, hjust = 0.5)
  )

# PD+ group
p_model_simulation_params_pd_on <- 
  ggplot(sim_means, aes(x = beta_log10, y = tau_log10, fill = meanReward)) +
  geom_raster() +
  scale_x_continuous(breaks = log10(bx), labels = bx, expand = c(0,0)) +
  scale_y_continuous(breaks = log10(by), labels = by, expand = c(0,0)) +
  labs(x = expression(paste('Exploration bonus ', beta)),
       y = expression(paste('Random exploration ', tau))) +
  # scale_fill_gradientn(colours = colorspace::sequential_hcl(500, "plasma"),name = "Normalized\nreward") +
  # scale_fill_viridis_c(option = "plasma", name = "Normalized\nreward") +
   scale_fill_gradientn(colours = cet_pal(5, name = "l7"), name = "Normalized\nreward") + #l3 l16 i5 "rainbow"
  coord_cartesian(xlim = c(-3,0.5), ylim=c(-3, 0.5)) +
  geom_jitter( # add individual points
    data = subset(marker2, group == "PD+"),
    aes(x = beta_log10, y = tau_log10, shape = group, colour = group),
    inherit.aes = FALSE,  
    size = 3,
    color = "white",
    fill = "#1b9e77",
    shape = 24,
    stroke = 0.3,
    width = 0.1
  ) +
 geom_point( # add median
    data = subset(marker, group == "PD+"),
    aes(x = beta_log10, y = tau_log10, colour = group),
    inherit.aes = FALSE,
    shape = 24,        
    size = 5,          
    stroke = 1.5,
    fill = "#1b9e77",
    color = "white"
  ) +
  ggtitle("PD+") +
  theme_classic() +
  theme(strip.background = element_blank(),  
        strip.text = element_text(color = "black", size=18),
        legend.position = "none",
        #        legend.justification = c(0, 1),
        axis.title = element_text(size = 18),
        axis.text = element_text(size = 18),
        plot.title = element_text(size = 18, hjust = 0.5),
        axis.title.y = element_blank()
        # axis.text.y = element_blank()
  )
# p_model_simulation_params_pd_on

# Extract legend
shared_legend <- cowplot::get_legend(
  p_model_simulation_params_control +
    theme(legend.position = "right",
          legend.title = element_text(size = 16),
          legend.text  = element_text(size = 14))
)

# combine plots
p_model_simulation_params_combined <- cowplot::plot_grid(
  p_model_simulation_params_pd_off,
  p_model_simulation_params_pd_on,
  p_model_simulation_params_control,
  nrow = 1, align = "hv", axis = "tblr", rel_widths = c(1,1,1),
  labels = NULL
)

# add  legend
p_model_simulation_params_combined <- cowplot::plot_grid(
  p_model_simulation_params_combined, shared_legend,
  ncol = 2, rel_widths = c(1, 0.10)   
)

# add title 
p_model_simulation_params <- ggdraw() +
  draw_label(
    "Model performance simulation",   
    x = 0, y = 0.98,               
    hjust = 0.5, vjust = 1,
    size = 24
  ) +
  draw_plot(p_model_simulation_params_combined, y = 0, height = 0.9)  

# p_model_simulation_params

ggsave("plots/model_simulation_params.png", p_model_simulation_params, width = 14, height = 5, dpi = 300)



# zoomed in version
# combine plots
p_model_simulation_params_combined_zoom <- cowplot::plot_grid(
  p_model_simulation_params_pd_off + coord_cartesian(xlim = c(-2,0.1), ylim=c(-2, -0.5)), 
  p_model_simulation_params_pd_on + coord_cartesian(xlim = c(-2,0.1), ylim=c(-2, -0.5)), 
  p_model_simulation_params_control + coord_cartesian(xlim = c(-2,0.1), ylim=c(-2, -0.5)) ,
  nrow = 1, align = "hv", axis = "tblr", rel_widths = c(1,1,1),
  labels = NULL
)

# add  legend
p_model_simulation_params_combined_zoom <- cowplot::plot_grid(
  p_model_simulation_params_combined_zoom, shared_legend,
  ncol = 2, rel_widths = c(1, 0.10)   
)

# add title 
p_model_simulation_params_zoom <- ggdraw() +
  draw_label(
    "Model performance simulation",   
    x = 0, y = 0.98,               
    hjust = 0, vjust = 1,
    size = 24
  ) +
  draw_plot(p_model_simulation_params_combined_zoom, y = 0, height = 0.9)  

p_model_simulation_params_zoom
ggsave("plots/model_simulation_params_zoom.png", p_model_simulation_params_zoom, width = 14, height = 5, dpi = 300)

Figure 4: Model simulation.

7 Supplementary Information

7.1 Statistical analyses

Statistical analyses were performed using R. We report both frequentist and Bayesian statistics, using Bayes factors (BF) to quantify the relative evidence of the data in favor of the alternative hypothesis (\(H_1\)) over the null (\(H_0\)). All data and code required for reproducing the statistical analyses and figures are available at ADD GITHUB or OSF LINK.

For parametric group comparisons, we report (paired or independent) Student’s t-tests (two-tailed). For non-parametric comparisons we used the Mann-Whitney U test or Wilcoxon signed-rank test. Bayes factors for the t-tests were computed with the package (Morey & Rouder, 2024), using its default settings. Bayes factor for rank tests were computed following (Doorn et al., 2020).

Linear correlations were assessed using Pearson’s \(r\), with the Bayes factors computed with the BayesFactor package (Morey & Rouder, 2024), using its default settings. Bayes factors for rank correlations quantified with Kendall’s tau were computed using an implementation from Doorn et al. (2018).

7.2 Supplementary computational results

8 Article figure

The following code generates Figure 3 from the article.

8.1 Figure 3. Computational results

# Combine
cowplot::plot_grid(
  p_model_comparison,
  p_classification_participants,
  p_gpucb_params,
  # p_parameters_reward,
  p_model_simulation_params_zoom,
  ncol = 1,
  # rel_heights = c(1.2, 1, 1),
  labels = c("auto"),
  label_y = 1.02, 
  label_size = 22,
  align = "v",     
  axis = "l"       
)

ggsave("plots/computational_results.png", dpi=300, height=16, width=14)

# # Figure inclduing GP-UCB visulaiztaion
# # get GP-UCB model illustration (Giron et al., 2023, NHB)
# img <- image_read("img/GP-UCB model.png")  
# gimg <- rasterGrob(as.raster(img), interpolate = TRUE)
# 
# # add title
# gp_ucp_model <- ggdraw() +
#   draw_label("Gaussian Process Upper Confidence Bound (GP-UCB) Model", size = 24, x = 0.05, y = 1, hjust = 0, vjust = 1.5) +
#   # draw_grob(gimg, x = 0.5, y = 0.5, width = 1, height = 1)
#   draw_grob(gimg)
# 
# # Combine
# cowplot::plot_grid(
#   gp_ucp_model,
#   p_model_comparison,
#   p_gpucb_params,             
#   ncol = 1,
#   # rel_heights = c(1.2, 1, 1),
#   labels = c("AUTO"),
#   label_size = 22
# )
# 
# ggsave("plots/computational_results.png", dpi=300, height=16, width=15)

Figure 5: Computational results.

9 Relations of model parameters and clinical indicators

9.1 Bivariate correlations

# 
df_params_clinical_indicators <- df_gpucb_params %>% 
  left_join(dat_sample, by= c("id", "group")) %>% 
  select(id, group, param, estimate, hoehn_yahr, BDI, MMSE) %>%
  pivot_wider(names_from = param, values_from = estimate) %>% 
  ungroup()


# control group
cor_matrix_controls <- df_params_clinical_indicators %>% 
  filter(group == "Control") %>% 
  select(BDI, MMSE, lambda, beta, tau) %>% 
  cor(method = "kendall", use = "pairwise.complete.obs")

p_cor_matrix_controls <- 
  ggcorrplot(cor_matrix_controls,
           # method = "circle",     
           type   = "lower",      
           lab    = TRUE,         
           show.legend = FALSE,    
           title  = "Controls") +
    scale_fill_gradient2(
    low = "darkred",   
    mid = "white",     
    high = "darkgreen",
    midpoint = 0,
    limit = c(-1, 1)
  ) +
  theme_classic() +
  theme(legend.position = "none",
        axis.text = element_text(size = 14),
        plot.title = element_text(size = 16, hjust = 0.5),
        axis.title = element_blank()
  ) +
  scale_x_discrete("",
    labels = c(
      "lambda" = "λ",
      "beta"   = "β",
      "tau"    = "τ",
      "hoehn_yahr" = "Hoehn–Yahr",
      "BDI"    = "BDI",
      "MMSE"   = "MMSE"
    )) +
  scale_y_discrete("",
    labels = c(
      "lambda" = "λ",
      "beta"   = "β",
      "tau"    = "τ",
      "hoehn_yahr" = "Hoehn–Yahr",
      "BDI"    = "BDI",
      "MMSE"   = "MMSE"
    ))

# PD+ patients on medication 
cor_matrix_pd_plus <- df_params_clinical_indicators %>% 
  filter(group == "PD+") %>% 
  select(hoehn_yahr, BDI, MMSE, lambda, beta, tau) %>% 
  cor(method = "pearson", use = "pairwise.complete.obs")

p_cor_matrix_pd_plus <- 
  ggcorrplot(cor_matrix_pd_plus,
           # method = "circle",     
           type   = "lower",      
           lab    = TRUE,         
           show.legend = FALSE,    
           title  = "PD+") +
    scale_fill_gradient2(
    low = "darkred",   
    mid = "white",     
    high = "darkgreen",
    midpoint = 0,
    limit = c(-1, 1)
  ) +
  theme_classic() +
  theme(legend.position = "none",
        axis.text = element_text(size = 14),
        plot.title = element_text(size = 16, hjust = 0.5),
        axis.title = element_blank()
  ) +
  scale_x_discrete("",
    labels = c(
      "lambda" = "λ",
      "beta"   = "β",
      "tau"    = "τ",
      "hoehn_yahr" = "Hoehn–Yahr",
      "BDI"    = "BDI",
      "MMSE"   = "MMSE"
    )) +
  scale_y_discrete("",
    labels = c(
      "lambda" = "λ",
      "beta"   = "β",
      "tau"    = "τ",
      "hoehn_yahr" = "Hoehn–Yahr",
      "BDI"    = "BDI",
      "MMSE"   = "MMSE"
    ))



# PD- patients off medication 
cor_matrix_pd_minus <- df_params_clinical_indicators %>% 
  filter(group == "PD-") %>% 
  select(hoehn_yahr, BDI, MMSE, lambda, beta, tau) %>% 
  cor(method = "pearson", use = "pairwise.complete.obs")

p_cor_matrix_pd_minus <- 
  ggcorrplot(cor_matrix_pd_minus,
           # method = "circle",     
           type   = "lower",      
           lab    = TRUE,         
           show.legend = FALSE,    
           title  = "PD-") +
    scale_fill_gradient2(
    low = "darkred",   
    mid = "white",     
    high = "darkgreen",
    midpoint = 0,
    limit = c(-1, 1)
  ) +
  theme_classic() +
  theme(legend.position = "none",
        axis.text = element_text(size = 14),
        plot.title = element_text(size = 16, hjust = 0.5),
        axis.title = element_blank()
  ) +
  scale_x_discrete("",
    labels = c(
      "lambda" = "λ",
      "beta"   = "β",
      "tau"    = "τ",
      "hoehn_yahr" = "Hoehn–Yahr",
      "BDI"    = "BDI",
      "MMSE"   = "MMSE"
    )) +
  scale_y_discrete("",
    labels = c(
      "lambda" = "λ",
      "beta"   = "β",
      "tau"    = "τ",
      "hoehn_yahr" = "Hoehn–Yahr",
      "BDI"    = "BDI",
      "MMSE"   = "MMSE"
    ))

p_cor_matrix <- cowplot::plot_grid(
  p_cor_matrix_pd_minus,
  p_cor_matrix_pd_plus,
  p_cor_matrix_controls,
  nrow = 1, align = "hv", axis = "tblr", rel_widths = c(1,1,1),
  labels = NULL
)

p_cor_matrix2 <- p_cor_matrix + draw_figure_label(label = "Correlation of model parameters and clinical indicators (Kendall's τ)", position = "top.left", size = 18)

ggsave("plots/cor_matrix_params_indicators.png", p_cor_matrix2, width = 12, height = 4)

9.2 Regression analyses

GP-UCB model parameters \(\lambda\) (amount of generalization), \(\beta\) (exploration bonus), and \(\tau\) (amount of random exploration) acorss allparticipants

# for now, random intercepts only, Random intercept + random slope not stable

# fit models: main effects only
lm_tau_bdi_mmse <- lm(log(tau) ~ group + BDI + MMSE , 
                             data = df_params_clinical_indicators)

lm_beta_bdi_mmse <- lm(log(beta) ~ group + BDI + MMSE , 
                             data = df_params_clinical_indicators)

lm_lambda_bdi_mmse <- lm(lambda ~ group + BDI + MMSE , 
                             data = df_params_clinical_indicators)


tab_model(lm_lambda_bdi_mmse, lm_beta_bdi_mmse, lm_tau_bdi_mmse)

	lambda			log(beta)			log(tau)
Predictors	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p
(Intercept)	-0.52	-2.14 – 1.10	0.526	10.81	-1.77 – 23.40	0.091	0.12	-14.16 – 14.40	0.987
group [PD+]	0.01	-0.10 – 0.12	0.893	-1.00	-1.85 – -0.16	0.021	-0.67	-1.63 – 0.29	0.169
group [Control]	0.10	-0.01 – 0.21	0.068	-1.41	-2.24 – -0.57	0.001	-0.62	-1.57 – 0.33	0.199
BDI	0.00	-0.01 – 0.02	0.624	-0.04	-0.13 – 0.06	0.452	-0.01	-0.12 – 0.10	0.817
MMSE	0.04	-0.02 – 0.09	0.198	-0.35	-0.78 – 0.07	0.104	-0.09	-0.57 – 0.40	0.724
Observations	97			97			97
R2 / R2 adjusted	0.060 / 0.020			0.146 / 0.109			0.029 / -0.013

# fit models: main effects + interactions with group
lm_tau_bdi_mmse <- lm(log(tau) ~ group *(BDI + MMSE), 
                             data = df_params_clinical_indicators)

lm_beta_bdi_mmse <- lm(log(beta) ~ group * (BDI + MMSE), 
                             data = df_params_clinical_indicators)

lm_lambda_bdi_mmse <- lm(lambda ~ group * (BDI + MMSE), 
                             data = df_params_clinical_indicators)


tab_model(lm_lambda_bdi_mmse, lm_beta_bdi_mmse, lm_tau_bdi_mmse)

	lambda			log(beta)			log(tau)
Predictors	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p
(Intercept)	0.24	-2.85 – 3.32	0.880	3.34	-19.40 – 26.08	0.771	-7.03	-33.87 – 19.82	0.604
group [PD+]	-1.38	-5.74 – 2.99	0.532	2.14	-30.00 – 34.27	0.895	2.68	-35.26 – 40.62	0.889
group [Control]	-0.82	-4.87 – 3.23	0.689	20.50	-9.37 – 50.37	0.176	17.19	-18.08 – 52.45	0.335
BDI	0.00	-0.02 – 0.02	0.925	-0.12	-0.29 – 0.04	0.135	0.00	-0.19 – 0.20	0.973
MMSE	0.01	-0.09 – 0.11	0.844	-0.07	-0.84 – 0.70	0.858	0.16	-0.75 – 1.07	0.732
group [PD+] × BDI	0.01	-0.03 – 0.04	0.726	0.07	-0.17 – 0.30	0.560	-0.07	-0.35 – 0.20	0.598
group [Control] × BDI	-0.00	-0.03 – 0.03	0.984	0.26	0.02 – 0.49	0.031	0.07	-0.21 – 0.34	0.641
group [PD+] × MMSE	0.05	-0.10 – 0.19	0.534	-0.13	-1.21 – 0.95	0.811	-0.10	-1.37 – 1.18	0.881
group [Control] × MMSE	0.03	-0.11 – 0.17	0.646	-0.83	-1.85 – 0.18	0.107	-0.64	-1.84 – 0.57	0.296
Observations	97			97			97
R2 / R2 adjusted	0.066 / -0.019			0.234 / 0.165			0.057 / -0.029

Linear regression with model GP-UCB model paramteres as function of group and depression level (BDI-II) and cognitive functioning (MMSE). All participants

GP-UCB model parameters \(\lambda\) (amount of generalization), \(\beta\) (exploration bonus), and \(\tau\) (amount of random exploration); PD patients only.

# for now, random intercepts only, Random intercept + random slope not stable

# fit models: main effects only
lm_tau_bdi_mmse_hy <- lm(log(tau) ~ group + BDI + MMSE + hoehn_yahr,
                             data = subset(df_params_clinical_indicators, group != "Control"))

lm_beta_bdi_mmse_hy <- lm(log(beta) ~ group + BDI + MMSE + hoehn_yahr,
                             data = subset(df_params_clinical_indicators, group != "Control"))

lm_lambda_bdi_mmse_hy <- lm(lambda ~ group + BDI + MMSE + hoehn_yahr, 
                             data = subset(df_params_clinical_indicators, group != "Control"))

tab_model(lm_lambda_bdi_mmse_hy, lm_beta_bdi_mmse_hy, lm_tau_bdi_mmse_hy)

	lambda			log(beta)			log(tau)
Predictors	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p
(Intercept)	-0.74	-2.88 – 1.39	0.487	7.04	-11.08 – 25.15	0.440	-5.30	-26.64 – 16.04	0.621
group [PD+]	0.01	-0.10 – 0.12	0.853	-1.08	-1.99 – -0.17	0.021	-0.73	-1.80 – 0.35	0.180
BDI	0.00	-0.01 – 0.02	0.688	-0.08	-0.21 – 0.04	0.195	-0.03	-0.18 – 0.12	0.676
MMSE	0.04	-0.03 – 0.11	0.256	-0.19	-0.79 – 0.42	0.542	0.11	-0.61 – 0.82	0.765
hoehn yahr	0.05	-0.04 – 0.13	0.282	-0.35	-1.07 – 0.36	0.329	0.01	-0.83 – 0.85	0.980
Observations	64			64			64
R2 / R2 adjusted	0.037 / -0.028			0.137 / 0.079			0.036 / -0.029

# fit models: main effects and interactions
lm_tau_bdi_mmse_hy <- lm(log(tau) ~ group * (BDI + MMSE + hoehn_yahr),
                             data = subset(df_params_clinical_indicators, group != "Control"))

lm_beta_bdi_mmse_hy <- lm(log(beta) ~ group * (BDI + MMSE + hoehn_yahr),
                             data = subset(df_params_clinical_indicators, group != "Control"))

lm_lambda_bdi_mmse_hy <- lm(lambda ~ group * (BDI + MMSE + hoehn_yahr), 
                             data = subset(df_params_clinical_indicators, group != "Control"))

tab_model(lm_lambda_bdi_mmse_hy, lm_beta_bdi_mmse_hy, lm_tau_bdi_mmse_hy)

	lambda			log(beta)			log(tau)
Predictors	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p
(Intercept)	-0.69	-3.62 – 2.24	0.640	6.50	-19.80 – 32.80	0.622	-10.59	-41.20 – 20.03	0.491
group [PD+]	-0.11	-4.25 – 4.04	0.959	0.07	-37.12 – 37.25	0.997	9.70	-33.59 – 52.99	0.655
BDI	-0.00	-0.02 – 0.02	0.782	-0.11	-0.29 – 0.07	0.237	-0.01	-0.23 – 0.20	0.913
MMSE	0.03	-0.06 – 0.13	0.501	-0.15	-1.02 – 0.73	0.738	0.24	-0.78 – 1.27	0.633
hoehn yahr	0.15	0.04 – 0.27	0.011	-0.53	-1.58 – 0.52	0.320	0.59	-0.63 – 1.82	0.335
group [PD+] × BDI	0.01	-0.02 – 0.04	0.494	0.06	-0.21 – 0.32	0.664	-0.05	-0.36 – 0.26	0.734
group [PD+] × MMSE	0.01	-0.12 – 0.15	0.829	-0.08	-1.32 – 1.15	0.897	-0.27	-1.71 – 1.17	0.709
group [PD+] × hoehn yahr	-0.21	-0.37 – -0.05	0.013	0.36	-1.11 – 1.82	0.628	-1.14	-2.84 – 0.56	0.185
Observations	64			64			64
R2 / R2 adjusted	0.145 / 0.039			0.148 / 0.041			0.070 / -0.046

Linear regression with model GP-UCB model paramteres as function of group and depression level (BDI-II) and cognitive functioning (MMSE). PD patients only.

10 Model params of participants best explained by GP-UCB model

For this analysis we only consider participants who were best explained by the GP-UCB model. The results are consistent with the same analyses performed with the full sample above: no substantial differences in amount of generalization \(\lambda\), marked differences in terms of the exploration bonus \(\beta\), and no differences in terms of random exploration \(\tau\). The only difference is that we found a difference between the control and off-medication group in the extent of generalization when using the full sample, whereas we found no difference when only considering the subset of participants best accounted for by the GP-UCB model.

10.0.1 Generalization \(\lambda\)

The parameter \(\lambda\) represents the length-scale in the RBF kernel, which governs the amount of generalization, i.e., to what extent participants assume a spatial correlation between options (higher \(\lambda\) = stronger generalization). Overall, the amount of generalization was very similar between groups.

• Control vs. PD+: \(U=210\), \(p=.402\), \(r_{ au}=.12\), \(BF=.44\)
• Control vs. PD-: \(U=219\), \(p=.150\), \(r_{ au}=.20\), \(BF=.63\)
• PD+ vs. PD-: \(U=191\), \(p=.558\), \(r_{ au}=.08\), \(BF=.34\)

10.0.2 Exploration bonus \(\beta\)

The parameter \(\beta\) represents the uncertainty bonus, i.e. how much expected rewards are positively inflated by their uncertainty (higher \(\beta\) = more uncertainty-directed exploration). Controls and PD+ patients on medication did not differ, and both groups had lower beta estimates than the dopamine-depleted patients in the PD− group. These differences suggest that levodopa medication modulated the amount of uncertainty-directed exploration by restoring beta to levels comparable to those observed in controls without PD. This aligns with findings from a restless bandit paradigm, where L-Dopa reduced the amount of directed exploration in healthy volunteers, while the level of random exploration remained unaffected (Chakroun et al., 2020).

• Control vs. PD+: \(U=182\), \(p=.977\), \(r_{ au}=.01\), \(BF=.32\)
• Control vs. PD-: \(U=55\), \(p<.001\), \(r_{ au}=-.49\), \(BF=14\)
• PD+ vs. PD-: \(U=54\), \(p<.001\), \(r_{ au}=-.49\), \(BF=31\)

10.0.3 Random exploration \(\tau\)

The parameter \(\tau\) represents the amount of decision noise, i.e. stochastic variability in the softmax decision rule (lower \(\tau\) = more decision noise, i.e. more uniform distribution; conversely, \(\tau \rightarrow \infty \quad \Rightarrow \quad \text{argmax (greedy)}\)). There were no group differences in rge temperature paramter \(\tau\), indicating comparable amounts of random exploration regardless of group.

• Control vs. PD+: \(U=193\), \(p=.729\), \(r_{ au}=.05\), \(BF=.35\)
• Control vs. PD-: \(U=162\), \(p=.799\), \(r_{ au}=-.04\), \(BF=.33\)
• PD+ vs. PD-: \(U=140\), \(p=.358\), \(r_{ au}=-.13\), \(BF=.44\)

df_gpucb_params_subset <- modelFits %>% 
  filter(best_ModelName == "GP-UCB") %>% 
  filter(kernel=='GP' & acq == 'UCB') %>% 
  pivot_longer(c('lambda', 'beta', 'tau'), names_to = 'param', values_to = 'estimate') %>% 
  mutate(param = factor(param, levels = c('lambda', 'beta', 'tau'))) 


ggboxplot(df_gpucb_params_subset, 
          x = "group", 
          y = "estimate",
          color = "group", palette =groupcolors, fill = "group", alpha = 0.2,
          add = "jitter", jitter.size = 0.5, shape = "group", title = "GP-UCB parameter estimates") +
  facet_wrap(~param, nrow = 1) +
  scale_y_log10(breaks = c(0.01, 0.1, 1, 10, 100), labels = c("0.01", "0.1", "1", "10", "100")) +
  ylab("Estimate (log scale") +
  xlab("") +
  stat_summary(fun = mean, geom="point", shape = 23, fill = "white", size=2) +
  theme_classic() +
  theme(strip.background = element_blank(),  
        strip.text = element_text(color = "black", size=12),
        legend.position = "none",
        plot.title = element_text(size = 18)
  )
  
# ggsave("plots/GP-UCB_params_subset.png", width = 9, height = 5)

Figure 6: Parameter estimates of GP-UCB model, estimated through leave-one-round-out cross validation. Each dot is one participant. Only participants are included who were best described by the GP-UCB model.

11 Session Information

sessionInfo()

R version 4.5.0 (2025-04-11)
Platform: aarch64-apple-darwin20
Running under: macOS Sequoia 15.6.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Europe/Berlin
tzcode source: internal

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] ggcorrplot_0.1.4.1     cetcolor_0.2.0         magick_2.8.7          
 [4] rstatix_0.7.2          parameters_0.27.0      ggthemes_5.1.0        
 [7] waffle_1.0.2           cowplot_1.2.0          scales_1.4.0          
[10] hms_1.1.3              ggpubr_0.6.1           viridis_0.6.5         
[13] viridisLite_0.4.2      emmeans_1.11.2         afex_1.4-1            
[16] kableExtra_1.4.0       brms_2.22.0            Rcpp_1.1.0            
[19] lsr_0.5.2              sjPlot_2.9.0           lme4_1.1-37           
[22] RColorBrewer_1.1-3     lubridate_1.9.4        forcats_1.0.0         
[25] stringr_1.5.1          dplyr_1.1.4            purrr_1.1.0           
[28] readr_2.1.5            tidyr_1.3.1            tibble_3.3.0          
[31] ggplot2_3.5.2          tidyverse_2.0.0        BayesFactor_0.9.12-4.7
[34] Matrix_1.7-3           coda_0.19-4.1          gridExtra_2.3         

loaded via a namespace (and not attached):
 [1] Rdpack_2.6.4         pbapply_1.7-4        rlang_1.1.6         
 [4] magrittr_2.0.3       snakecase_0.11.1     matrixStats_1.5.0   
 [7] compiler_4.5.0       mgcv_1.9-3           loo_2.8.0           
[10] systemfonts_1.2.3    vctrs_0.6.5          reshape2_1.4.4      
[13] crayon_1.5.3         pkgconfig_2.0.3      fastmap_1.2.0       
[16] backports_1.5.0      labeling_0.4.3       effectsize_1.0.1    
[19] rmarkdown_2.29       tzdb_0.5.0           nloptr_2.2.1        
[22] ragg_1.4.0           bit_4.6.0            MatrixModels_0.5-4  
[25] xfun_0.52            jsonlite_2.0.0       sjmisc_2.8.11       
[28] broom_1.0.9          parallel_4.5.0       R6_2.6.1            
[31] stringi_1.8.7        car_3.1-3            boot_1.3-31         
[34] extrafontdb_1.0      numDeriv_2016.8-1.1  estimability_1.5.1  
[37] knitr_1.50           extrafont_0.19       bayesplot_1.13.0    
[40] splines_4.5.0        timechange_0.3.0     tidyselect_1.2.1    
[43] rstudioapi_0.17.1    abind_1.4-8          yaml_2.3.10         
[46] sjlabelled_1.2.0     curl_6.4.0           lattice_0.22-7      
[49] lmerTest_3.1-3       plyr_1.8.9           bayestestR_0.16.1   
[52] withr_3.0.2          bridgesampling_1.1-2 posterior_1.6.1     
[55] evaluate_1.0.4       RcppParallel_5.1.10  xml2_1.3.8          
[58] pillar_1.11.0        carData_3.0-5        tensorA_0.36.2.1    
[61] checkmate_2.3.2      DT_0.33              insight_1.3.1       
[64] reformulas_0.4.1     distributional_0.5.0 generics_0.1.4      
[67] vroom_1.6.5          rstantools_2.4.0     minqa_1.2.8         
[70] glue_1.8.0           tools_4.5.0          ggsignif_0.6.4      
[73] mvtnorm_1.3-3        Rttf2pt1_1.3.12      rbibutils_2.3       
[76] datawizard_1.2.0     nlme_3.1-168         performance_0.15.0  
[79] Formula_1.2-5        cli_3.6.5            textshaping_1.0.1   
[82] svglite_2.2.1        Brobdingnag_1.2-9    gtable_0.3.6        
[85] digest_0.6.37        htmlwidgets_1.6.4    farver_2.1.2        
[88] htmltools_0.5.8.1    lifecycle_1.0.4      bit64_4.6.0-1       
[91] MASS_7.3-65         

References

Chakroun, K., Mathar, D., Wiehler, A., Ganzer, F., & Peters, J. (2020). Dopaminergic modulation of the exploration/exploitation trade-off in human decision-making. eLife, 9, e51260.

Doorn, J. van, Ly, A., Marsman, M., & Wagenmakers, E.-J. (2018). Bayesian inference for kendall’s rank correlation coefficient. The American Statistician, 72, 303–308.

Doorn, J. van, Ly, A., Marsman, M., & Wagenmakers, E.-J. (2020). Bayesian rank-based hypothesis testing for the rank sum test, the signed rank test, and spearman’s \(\rho\). Journal of Applied Statistics, 47(16), 2984–3006.

Giron, A. P., Ciranka, S., Schulz, E., Bos, W. van den, Ruggeri, A., Meder, B., & Wu, C. M. (2023). Developmental changes in exploration resemble stochastic optimization. Nature Human Behaviour, 7(11), 1955–1967. https://doi.org/https://doi.org/10.1038/s41562-023-01662-1

Meder, B., Wu, C. M., Schulz, E., & Ruggeri, A. (2021). Development of directed and random exploration in children. Developmental Science, 24(4), e13095. https://doi.org/https://doi.org/10.1111/desc.13095

Morey, R. D., & Rouder, J. N. (2024). BayesFactor: Computation of bayes factors for common designs. https://doi.org/10.32614/CRAN.package.BayesFactor

Schulz, E., Wu, C. M., Ruggeri, A., & Meder, B. (2019). Searching for rewards like a child means less generalization and more directed exploration. Psychological Science, 30(11), 1561–1572. https://doi.org/10.1177/0956797619863663

Wu, C. M., Meder, B., & Schulz, E. (2025). Unifying principles of generalization: Past, present, and future. Annual Review of Psychology, 76, 275–302. https://doi.org/https://doi.org/10.1146/annurev-psych-021524-110810

Wu, C. M., Schulz, E., Garvert, M. M., Meder, B., & Schuck, N. W. (2020). Similarities and differences in spatial and non-spatial cognitive maps. PLOS Computational Biology, 16(9), e1008149. https://doi.org/10.1371/journal.pcbi.1008149

Wu, C. M., Schulz, E., Speekenbrink, M., Nelson, J. D., & Meder, B. (2018). Generalization guides human exploration in vast decision spaces. Nature Human Behaviour, 2, 915–924. https://doi.org/10.1038/s41562-018-0467-4